HerbPt1 16S Figures

require(phyloseq)
require(tidyverse)
require(phyloseq)
require(reshape2)
require(dplyr)
require(ggplot2)
require(microbiome)

Loading required package: microbiome

microbiome R package (microbiome.github.com)
    


 Copyright (C) 2011-2020 Leo Lahti, 
    Sudarshan Shetty et al. <microbiome.github.io>


Attaching package: ‘microbiome’

The following object is masked from ‘package:micrUBIfuns’:

    prevalence

The following object is masked from ‘package:psych’:

    alpha

The following object is masked from ‘package:vegan’:

    diversity

The following object is masked from ‘package:ggplot2’:

    alpha

The following object is masked from ‘package:base’:

    transform

require(vegan)

Load data

ps_dmn <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/ITS/DMN_ests_ITS.Rdata")
sample_data(ps_dmn)$Herbicide <- factor(sample_data(ps_dmn)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_dmn)$herb_time<-paste(sample_data(ps_dmn)$Herbicide, sample_data(ps_dmn)$Time, sep = "_")

#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_dmn)$Mode<-sample_data(ps_dmn)$Herbicide

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")

sample_data(ps_dmn)$Mode<- as.factor(values[match(sample_data(ps_dmn)$Mode, index)])

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")

sample_data(ps_dmn)$Herbicide <- as.factor(values[match(sample_data(ps_dmn)$Herbicide, index)])


ps_rare <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/ITS/HerbPt1_rare_ITS.Rdata")
sample_data(ps_rare)$Herbicide <- factor(sample_data(ps_rare)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_rare)$herb_time<-paste(sample_data(ps_rare)$Herbicide, sample_data(ps_rare)$Time, sep = "_")


#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")

sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")

sample_data(ps_rare)$Herbicide <- as.factor(values[match(sample_data(ps_rare)$Herbicide, index)])


ps_trans <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/ITS/HerbPt1_hel_trans_ITS.Rdata")
sample_data(ps_trans)$Herbicide <- factor(sample_data(ps_trans)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_trans)$herb_time<-paste(sample_data(ps_trans)$Herbicide, sample_data(ps_trans)$Time, sep = "_")


#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_trans)$Mode<-sample_data(ps_trans)$Herbicide

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")

sample_data(ps_trans)$Mode<- as.factor(values[match(sample_data(ps_trans)$Mode, index)])

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")

sample_data(ps_trans)$Herbicide <- as.factor(values[match(sample_data(ps_trans)$Herbicide, index)])

Remove samples that are outlines or under sequenced.

ps_dmn <-  subset_samples(ps_dmn, sample_names(ps_dmn) != "G009SG")
ps_dmn <-  subset_samples(ps_dmn, sample_names(ps_dmn) != "G095SG")
ps_dmn <-  subset_samples(ps_dmn, sample_names(ps_dmn) != "G123SG")
ps_dmn <-  subset_samples(ps_dmn, sample_names(ps_dmn) != "G129SG")

ps_rare <-  subset_samples(ps_rare, sample_names(ps_rare) != "G009SG")
ps_rare <-  subset_samples(ps_rare, sample_names(ps_rare) != "G095SG")
ps_rare <-  subset_samples(ps_rare, sample_names(ps_rare) != "G123SG")
ps_rare <-  subset_samples(ps_rare, sample_names(ps_rare) != "G129SG")
ps_rare_sub<-prune_taxa(taxa_sums(ps_rare) > 2, ps_rare)

ordinations and adonis testing with three separate objects (i.e., dmn, rarefied, transformed). Rare taxa are removed from rarefied and transformed to successfully ordinate. At this point, the transformed data will not ordinate. This section is full dataset ordinations.

ord_dmn<-ordinate(physeq = ps_dmn, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit = 1000)

ord_rare<-ordinate(physeq = ps_rare_sub, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit = 1000)

ps_trans_sub<-prune_taxa(taxa_sums(ps_trans) > 0.01, ps_trans)
ord_transformed<-ordinate(physeq = ps_trans_sub, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit = 1000)

create alphadiversity tables

sample_sums(ps_rare)

G001SG G002SG G003SG G004SG G005SG G006SG G007SG G008SG G010SG G012SG G013SG G014SG G016SG G017SG G018SG G019SG G020SG 
  1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000 
G021SG G022SG G024SG G025SG G026SG G027SG G028SG G029SG G030SG G031SG G032SG G033SG G034SG G035SG G036SG G037SG G038SG 
  1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000 
G039SG G040SG G041SG G043SG G044SG G045SG G046SG G047SG G048SG G050SG G051SG G052SG G053SG G054SG G056SG G057SG G058SG 
  1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000 
G059SG G060SG G062SG G063SG G064SG G065SG G066SG G067SG G068SG G069SG G070SG G071SG G072SG G073SG G074SG G075SG G076SG 
  1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000 
G077SG G078SG G079SG G080SG G081SG G082SG G083SG G084SG G085SG G086SG G087SG G088SG G090SG G091SG G093SG G094SG G096SG 
  1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000 
G097SG G098SG G100SG G101SG G102SG G103SG G104SG G105SG G106SG G107SG G108SG G109SG G111SG G112SG G113SG G114SG G115SG 
  1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000 
G117SG G118SG G121SG G122SG G124SG G125SG G126SG G127SG G128SG G130SG G131SG G132SG G133SG G134SG G135SG G136SG G137SG 
  1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000 
G139SG G140SG G141SG G142SG G143SG G144SG G145SG G146SG G147SG G148SG G149SG G150SG G151SG G153SG G154SG G155SG G156SG 
  1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000 
G157SG G158SG G159SG G160SG G161SG G162SG G163SG G164SG G166SG G167SG G168SG G169SG G170SG G171SG G172SG G173SG G174SG 
  1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000   1000 
G176SG G178SG G179SG G180SG 
  1000   1000   1000   1000 

ps_rare

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 2033 taxa and 157 samples ]
sample_data() Sample Data:       [ 157 samples by 53 sample variables ]
tax_table()   Taxonomy Table:    [ 2033 taxa by 7 taxonomic ranks ]

alpha_div <- estimate_richness(physeq = ps_rare, measures = c("Observed", "Shannon", "Chao1"))
#pull out metadata and concatonate with alpha diversity metrics
md<-data.frame(sample_data(ps_rare))
alpha_div_md <- rownames_to_column(alpha_div, "Barcode_ID_G") %>% full_join(md) 

Joining, by = "Barcode_ID_G"

alpha_div_md$Herbicide <- factor(alpha_div_md$Herbicide, levels = c("Non-Treated", "Handweeded", "Atrazine-Mesotrione", "Dicamba", "Glyphosate"))

Shannon Div plots - no significant differences among herbicide treatments at any of the three time points

ggplot(data = alpha_div_md, aes(Herbicide, Shannon, color= Herbicide)) + facet_grid(. ~ Time) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 45, hjust = 1) )

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_Shannon.pdf")

Saving 7.29 x 4.51 in image

aov_t1<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T1",])
plot(aov_t1$residuals)

summary(aov_t1)

            Df Sum Sq Mean Sq F value Pr(>F)
Herbicide    4  0.389 0.09720    1.26  0.299
Residuals   48  3.704 0.07716               

aov_t2<-aov(Shannon~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T2",])
plot(aov_t2$residuals)

summary(aov_t2)

            Df Sum Sq Mean Sq F value Pr(>F)
Herbicide    4 0.1234 0.03086   0.508   0.73
Residuals   46 2.7936 0.06073               

aov_t3<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T3",])
plot(aov_t3$residuals)

summary(aov_t3)

            Df Sum Sq Mean Sq F value Pr(>F)  
Herbicide    4 0.5376 0.13440   2.551  0.051 .
Residuals   48 2.5287 0.05268                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(aov_t3, "Herbicide")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T3", ])

$Herbicide
                                        diff         lwr       upr     p adj
Handweeded-Non-Treated           0.086844630 -0.19738244 0.3710717 0.9078666
Atrazine-Mesotrione-Non-Treated  0.129708505 -0.14882205 0.4082391 0.6803851
Dicamba-Non-Treated             -0.002861275 -0.30174872 0.2960262 0.9999999
Glyphosate-Non-Treated           0.273363070 -0.01086400 0.5575901 0.0647237
Atrazine-Mesotrione-Handweeded   0.042863876 -0.22867317 0.3144009 0.9914461
Dicamba-Handweeded              -0.089705904 -0.38208716 0.2026754 0.9066040
Glyphosate-Handweeded            0.186518440 -0.09085878 0.4638957 0.3282566
Dicamba-Atrazine-Mesotrione     -0.132569780 -0.41941651 0.1542769 0.6864852
Glyphosate-Atrazine-Mesotrione   0.143654565 -0.12788248 0.4151916 0.5679081
Glyphosate-Dicamba               0.276224345 -0.01615691 0.5686056 0.0724032

Adonis testing of herbicide treatments by time point

ps_adonis<-function(physeq){
  otu_tab<-data.frame(phyloseq::otu_table(physeq))
  md_tab<-data.frame(phyloseq::sample_data(physeq))
    if(taxa_are_rows(physeq)== T){
       physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
            else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
  print(anova(vegan::betadisper(physeq_dist, md_tab$Herbicide)))
  vegan::adonis(physeq_dist ~ Herbicide * Time + Total_Weed_Veg , data = md_tab, permutations = 1000)
}

#ps_adonis(ps_rare_sub)
#remove one sample with no vegetation measurement. 
ps_rare_sub_57<-subset_samples(ps_rare_sub, sample_names(ps_rare_sub) != "G065SG")
ps_adonis(ps_rare_sub_57)
#ps_adonis(ps_trans_sub)
ps_dmn_57<-subset_samples(ps_dmn, sample_names(ps_dmn) != "G065SG")
#ps_adonis(ps_dmn)
ps_adonis(ps_dmn_57)

Ordination plots DMN by time point

ord_t1_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Run 0 stress 0.1335559 
Run 1 stress 0.1309245 
... New best solution
... Procrustes: rmse 0.06976384  max resid 0.2060325 
Run 2 stress 0.1315098 
Run 3 stress 0.132755 
Run 4 stress 0.134753 
Run 5 stress 0.1308758 
... New best solution
... Procrustes: rmse 0.0538814  max resid 0.2344595 
Run 6 stress 0.1331698 
Run 7 stress 0.1323904 
Run 8 stress 0.1346675 
Run 9 stress 0.1350022 
Run 10 stress 0.1322358 
Run 11 stress 0.1334733 
Run 12 stress 0.1315248 
Run 13 stress 0.1331407 
Run 14 stress 0.1327605 
Run 15 stress 0.1332143 
Run 16 stress 0.1369199 
Run 17 stress 0.1339829 
Run 18 stress 0.1331858 
Run 19 stress 0.1360529 
Run 20 stress 0.1352334 
Run 21 stress 0.1308759 
... Procrustes: rmse 0.001713885  max resid 0.006848176 
... Similar to previous best
*** Solution reached

T1_dmn<-ggordiplots::gg_ordiplot(ord = ord_t1_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T1_dmn$plot + theme_classic()

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_dmn_T1.pdf")

Saving 7.29 x 4.51 in image

ord_t2_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Run 0 stress 0.1325262 
Run 1 stress 0.1325269 
... Procrustes: rmse 0.0001074309  max resid 0.000408716 
... Similar to previous best
Run 2 stress 0.1325262 
... New best solution
... Procrustes: rmse 0.001015874  max resid 0.005488789 
... Similar to previous best
Run 3 stress 0.1327962 
... Procrustes: rmse 0.01038273  max resid 0.06195753 
Run 4 stress 0.1338926 
Run 5 stress 0.1392078 
Run 6 stress 0.1338929 
Run 7 stress 0.1325264 
... Procrustes: rmse 0.001066291  max resid 0.005401677 
... Similar to previous best
Run 8 stress 0.1325261 
... New best solution
... Procrustes: rmse 0.0001868934  max resid 0.001066974 
... Similar to previous best
Run 9 stress 0.1334238 
Run 10 stress 0.1325259 
... New best solution
... Procrustes: rmse 0.0007505591  max resid 0.004597761 
... Similar to previous best
Run 11 stress 0.132925 
... Procrustes: rmse 0.01264762  max resid 0.07728014 
Run 12 stress 0.1438152 
Run 13 stress 0.132526 
... Procrustes: rmse 0.0003451055  max resid 0.001743723 
... Similar to previous best
Run 14 stress 0.1325259 
... Procrustes: rmse 0.0003370863  max resid 0.00167214 
... Similar to previous best
Run 15 stress 0.1334866 
Run 16 stress 0.1325266 
... Procrustes: rmse 0.0009514394  max resid 0.005508805 
... Similar to previous best
Run 17 stress 0.1325265 
... Procrustes: rmse 0.0008955128  max resid 0.005161325 
... Similar to previous best
Run 18 stress 0.1449845 
Run 19 stress 0.1342762 
Run 20 stress 0.1325263 
... Procrustes: rmse 0.0003284059  max resid 0.001374251 
... Similar to previous best
*** Solution reached

T2_dmn<-ggordiplots::gg_ordiplot(ord = ord_t2_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T2_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_dmn_T2.pdf")

Saving 7.29 x 4.51 in image

#this time point needs to be checked out. The ordination is not working properly. 
ord_t3_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Run 0 stress 0.1451733 
Run 1 stress 0.1459287 
Run 2 stress 0.1510555 
Run 3 stress 0.1501906 
Run 4 stress 0.1503806 
Run 5 stress 0.1459215 
Run 6 stress 0.1528692 
Run 7 stress 0.1451733 
... Procrustes: rmse 0.00281782  max resid 0.01088482 
Run 8 stress 0.1502075 
Run 9 stress 0.1507344 
Run 10 stress 0.1478567 
Run 11 stress 0.1490906 
Run 12 stress 0.1451724 
... New best solution
... Procrustes: rmse 0.0002730257  max resid 0.001361369 
... Similar to previous best
Run 13 stress 0.1507808 
Run 14 stress 0.146477 
Run 15 stress 0.1510089 
Run 16 stress 0.1507665 
Run 17 stress 0.1518463 
Run 18 stress 0.1459425 
Run 19 stress 0.1451714 
... New best solution
... Procrustes: rmse 0.001928599  max resid 0.006972218 
... Similar to previous best
Run 20 stress 0.1459215 
*** Solution reached

T3_dmn<-ggordiplots::gg_ordiplot(ord = ord_t3_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T3_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_dmn_T3.pdf")

Saving 7.29 x 4.51 in image

Ordination plots on rarefied data by time point.

ord_t1_rare<-ordinate(physeq = subset_samples(ps_rare_sub, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.15407 
Run 1 stress 0.1540701 
... Procrustes: rmse 0.0004569647  max resid 0.001896947 
... Similar to previous best
Run 2 stress 0.1560276 
Run 3 stress 0.1540705 
... Procrustes: rmse 0.00348766  max resid 0.02182018 
Run 4 stress 0.1560277 
Run 5 stress 0.158209 
Run 6 stress 0.1586751 
Run 7 stress 0.1545097 
... Procrustes: rmse 0.01243307  max resid 0.04244993 
Run 8 stress 0.1540702 
... Procrustes: rmse 0.003479511  max resid 0.02173529 
Run 9 stress 0.1543042 
... Procrustes: rmse 0.01763231  max resid 0.08581444 
Run 10 stress 0.1541252 
... Procrustes: rmse 0.007485728  max resid 0.03683984 
Run 11 stress 0.1540706 
... Procrustes: rmse 0.003606795  max resid 0.02250138 
Run 12 stress 0.1548484 
Run 13 stress 0.1541275 
... Procrustes: rmse 0.00762305  max resid 0.03699455 
Run 14 stress 0.1543157 
... Procrustes: rmse 0.01699931  max resid 0.07997088 
Run 15 stress 0.1560286 
Run 16 stress 0.1560286 
Run 17 stress 0.1581542 
Run 18 stress 0.1540702 
... Procrustes: rmse 0.003594385  max resid 0.02243185 
Run 19 stress 0.1550319 
Run 20 stress 0.1560267 
*** Best solution repeated 1 times

T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T1_rare$plot + theme_classic()

T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T1_rare_plot<-T1_rare$plot + theme_classic() + xlim(-0.5, 0.5) + ylim(-0.5, 0.5)  + guides(color=guide_legend("Treatment")) + xlab("NMDS 1") + ylab("NMDS 2")
T1_rare_plot

Warning: Removed 3 rows containing missing values (`geom_point()`).

library(cowplot)
my_legend <- get_legend(T1_rare_plot)

Warning: Removed 3 rows containing missing values (`geom_point()`).

library(ggpubr)
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordinationlegend.pdf")

Saving 7.29 x 4.51 in image

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordinationlegend.eps")

Saving 7.29 x 4.51 in image

T1_rare_plot<-T1_rare_plot + theme(legend.position = "none")
T1_rare_plot

Warning: Removed 3 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T1.pdf")

Saving 7.29 x 4.51 in image

Warning: Removed 3 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T1.eps")

Saving 7.29 x 4.51 in image

Warning: Removed 3 rows containing missing values (`geom_point()`).

ord_t2_rare<-ordinate(physeq = subset_samples(ps_rare_sub, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1716481 
Run 1 stress 0.1720658 
... Procrustes: rmse 0.04180674  max resid 0.1226574 
Run 2 stress 0.1720724 
... Procrustes: rmse 0.0407305  max resid 0.1223211 
Run 3 stress 0.1701444 
... New best solution
... Procrustes: rmse 0.08690941  max resid 0.2978527 
Run 4 stress 0.1709893 
Run 5 stress 0.1711196 
Run 6 stress 0.1728237 
Run 7 stress 0.1708441 
Run 8 stress 0.1706909 
Run 9 stress 0.1706258 
... Procrustes: rmse 0.06561966  max resid 0.2792094 
Run 10 stress 0.171183 
Run 11 stress 0.1746478 
Run 12 stress 0.1746525 
Run 13 stress 0.1709914 
Run 14 stress 0.1734335 
Run 15 stress 0.1710644 
Run 16 stress 0.1710837 
Run 17 stress 0.1702885 
... Procrustes: rmse 0.05689984  max resid 0.2828329 
Run 18 stress 0.1734279 
Run 19 stress 0.1715313 
Run 20 stress 0.1713459 
Run 21 stress 0.1729146 
Run 22 stress 0.1736516 
Run 23 stress 0.1703011 
... Procrustes: rmse 0.05756458  max resid 0.2826533 
Run 24 stress 0.170112 
... New best solution
... Procrustes: rmse 0.005398775  max resid 0.02579526 
Run 25 stress 0.1702884 
... Procrustes: rmse 0.05862291  max resid 0.2821017 
Run 26 stress 0.1706282 
Run 27 stress 0.172058 
Run 28 stress 0.174658 
Run 29 stress 0.1706231 
Run 30 stress 0.1702879 
... Procrustes: rmse 0.05856164  max resid 0.2821153 
Run 31 stress 0.1701561 
... Procrustes: rmse 0.005990811  max resid 0.02698447 
Run 32 stress 0.1715783 
Run 33 stress 0.1757041 
Run 34 stress 0.1709007 
Run 35 stress 0.1715863 
Run 36 stress 0.1742448 
Run 37 stress 0.1702882 
... Procrustes: rmse 0.05868747  max resid 0.2820668 
Run 38 stress 0.1702889 
... Procrustes: rmse 0.05873721  max resid 0.2820447 
Run 39 stress 0.1715866 
Run 40 stress 0.170299 
... Procrustes: rmse 0.05905885  max resid 0.2819334 
Run 41 stress 0.1709006 
Run 42 stress 0.1702394 
... Procrustes: rmse 0.01169203  max resid 0.03765 
Run 43 stress 0.1767481 
Run 44 stress 0.1720567 
Run 45 stress 0.1711196 
Run 46 stress 0.1703014 
... Procrustes: rmse 0.05919363  max resid 0.2818825 
Run 47 stress 0.1709005 
Run 48 stress 0.1703017 
... Procrustes: rmse 0.05921868  max resid 0.2818621 
Run 49 stress 0.1709909 
Run 50 stress 0.1745794 
Run 51 stress 0.1709004 
Run 52 stress 0.1710641 
Run 53 stress 0.1731 
Run 54 stress 0.1728168 
Run 55 stress 0.1706225 
Run 56 stress 0.1706709 
Run 57 stress 0.1732789 
Run 58 stress 0.1711185 
Run 59 stress 0.1702091 
... Procrustes: rmse 0.009167803  max resid 0.02910266 
Run 60 stress 0.1703042 
... Procrustes: rmse 0.05940232  max resid 0.2817527 
Run 61 stress 0.1701124 
... Procrustes: rmse 0.0009526391  max resid 0.002988642 
... Similar to previous best
*** Best solution repeated 1 times

T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T2_rare_plot<-T2_rare$plot + theme_classic() + xlim(-0.5, 0.5) + ylim(-0.5, 0.5) + theme(legend.position = "none")  + xlab("NMDS 1") + ylab("NMDS 2")
T2_rare_plot

Warning: Removed 4 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T2.pdf")

Saving 7.29 x 4.51 in image

Warning: Removed 4 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T2.eps")

Saving 7.29 x 4.51 in image

Warning: Removed 4 rows containing missing values (`geom_point()`).

ord_t3_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1685562 
Run 1 stress 0.173731 
Run 2 stress 0.1770624 
Run 3 stress 0.1745126 
Run 4 stress 0.1722637 
Run 5 stress 0.1686696 
... Procrustes: rmse 0.01027886  max resid 0.05714972 
Run 6 stress 0.1722627 
Run 7 stress 0.1684508 
... New best solution
... Procrustes: rmse 0.01175199  max resid 0.06578898 
Run 8 stress 0.1786699 
Run 9 stress 0.1686779 
... Procrustes: rmse 0.02156585  max resid 0.130108 
Run 10 stress 0.1685888 
... Procrustes: rmse 0.01036983  max resid 0.05925135 
Run 11 stress 0.1769021 
Run 12 stress 0.1698586 
Run 13 stress 0.1798158 
Run 14 stress 0.1774968 
Run 15 stress 0.1723331 
Run 16 stress 0.1734184 
Run 17 stress 0.1704177 
Run 18 stress 0.1686749 
... Procrustes: rmse 0.02108582  max resid 0.1271145 
Run 19 stress 0.1686812 
... Procrustes: rmse 0.02201686  max resid 0.1329718 
Run 20 stress 0.1684506 
... New best solution
... Procrustes: rmse 0.0001109389  max resid 0.0005559503 
... Similar to previous best
*** Best solution repeated 1 times

T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T3_rare_plot<-T3_rare$plot + theme_classic() + xlim(-0.5, 0.5) + ylim(-0.5, 0.5) + theme(legend.position = "none")  + xlab("NMDS 1") + ylab("NMDS 2")
T3_rare_plot

Warning: Removed 4 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T3.pdf")

Saving 7.29 x 4.51 in image

Warning: Removed 4 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T3.eps")

Saving 7.29 x 4.51 in image

Warning: Removed 4 rows containing missing values (`geom_point()`).

library(ggpubr)
ggarrange(T1_rare_plot, T2_rare_plot, T3_rare_plot, ncol = 1)

Warning: Removed 3 rows containing missing values (`geom_point()`).
Warning: Removed 4 rows containing missing values (`geom_point()`).
Warning: Removed 4 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_ordination.pdf", width = 5, height = 10)

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_ordination.eps", width = 5, height = 10)

CAP ordination plots rarefied

t1_dist <- distance(subset_samples(ps_rare, Time=="T1"), method="bray") #get wUnifrac and save
t1_table<-as.matrix(dist(t1_dist)) #transform wUnifrac index
ord_t1_rare_cap <- capscale(t1_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T1"))))
T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T1_cap.pdf")


t2_dist <- distance(subset_samples(ps_rare, Time=="T2"), method="bray") #get wUnifrac and save
t2_table<-as.matrix(dist(t2_dist)) #transform wUnifrac index
ord_t2_rare_cap <- capscale(t2_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T2"))))
T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T2_cap.pdf")


#G166SG identified as outlier based on plots with it included. Removed to create plot. 
ps_rare <-  subset_samples(ps_rare, sample_names(ps_rare) != "G166SG")
t3_dist <- distance(subset_samples(ps_rare, Time=="T3"), method="bray") #get wUnifrac and save
t3_table<-as.matrix(dist(t3_dist)) #transform wUnifrac index
ord_t3_rare_cap <- capscale(t3_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T3"))))
T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T3_cap.pdf")

Pairwise adonis testing no needed becasue of insignificant gloabal test.

ps_pairwiseadonis<-function(physeq){
  otu_tab<-data.frame(phyloseq::otu_table(physeq))
  md_tab<-data.frame(phyloseq::sample_data(physeq))
    if(taxa_are_rows(physeq)== T){
       physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
            else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
pairwiseAdonis::pairwise.adonis(x = physeq_dist, factors = md_tab$Herbicide, p.adjust.m = "none", perm = 1000)
}

ps_t1<-subset_samples(ps_rare_sub, Time == "T1")
ps_t1<-prune_taxa(taxa_sums(ps_t1) > 2, ps_t1)

ps_t2<-subset_samples(ps_rare_sub, Time == "T2")
ps_t2<-prune_taxa(taxa_sums(ps_t2) > 2, ps_t2)

ps_t3<-subset_samples(ps_rare_sub, Time == "T3")
ps_t3<-prune_taxa(taxa_sums(ps_t3) > 2, ps_t3)


ps_pairwiseadonis(ps_t1)
ps_pairwiseadonis(ps_t2)
ps_pairwiseadonis(ps_t3)

Pairwise betadispr by treatment, time and mode

ps_betadispr<-function(physeq, groupingvar = "Groupingvar"){
  otu_tab<-data.frame(phyloseq::otu_table(physeq))
  md_tab<-data.frame(phyloseq::sample_data(physeq))
    if(taxa_are_rows(physeq)== T){
       physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
            else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
                mod<-vegan::betadisper(physeq_dist, md_tab[,groupingvar])
        ## Perform test
                print(anova(mod))
        ## Permutation test for F
                pmod <- vegan::permutest(mod, permutations = 1000, pairwise = TRUE)
                print(pmod)
                print(boxplot(mod))
}

permute test of dispersion

ps_betadispr(subset_samples(ps_rare_sub, Time == "T3"), groupingvar = "Herbicide")

Analysis of Variance Table

Response: Distances
          Df  Sum Sq   Mean Sq F value Pr(>F)
Groups     4 0.00635 0.0015867  0.1932 0.9408
Residuals 48 0.39424 0.0082133               

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     4 0.00635 0.0015867 0.1932   1000 0.9451
Residuals 48 0.39424 0.0082133                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                    Atrazine-Mesotrione Dicamba Glyphosate Handweeded Non-Treated
Atrazine-Mesotrione                     0.93007    0.53946    0.86913      0.7502
Dicamba                         0.92524            0.42557    0.73726      0.5854
Glyphosate                      0.54158 0.42422               0.62537      0.7832
Handweeded                      0.84582 0.73548    0.61029                 0.8442
Non-Treated                     0.71761 0.59620    0.77349    0.82944            
$stats
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.2201596 0.2491324 0.1949060 0.2643697 0.2208840
[2,] 0.2635563 0.2905369 0.2531917 0.2745103 0.2571909
[3,] 0.2896573 0.3043022 0.2848346 0.3109094 0.2923035
[4,] 0.3593706 0.3524222 0.3591617 0.3254448 0.3874820
[5,] 0.4966550 0.3757590 0.4464053 0.3722896 0.4658206

$n
[1] 12  9 11 11 10

$conf
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.2459558 0.2717092 0.2343518 0.2866449 0.2272048
[2,] 0.3333588 0.3368951 0.3353175 0.3351740 0.3574021

$out
   G141SG    G146SG    G144SG 
0.6092207 0.5282037 0.5427749 

$group
[1] 1 2 4

$names
[1] "Atrazine-Mesotrione" "Dicamba"             "Glyphosate"          "Handweeded"         
[5] "Non-Treated"        

beta_boxplot<-function (physeq, method = "bray", group) 
{
  require("phyloseq")
  require("ggplot2")
  group2samp <- list()
  group_list <- get_variable(sample_data(physeq), group)
  for (groups in levels(group_list)) {
    target_group <- which(group_list == groups)
    group2samp[[groups]] <- sample_names(physeq)[target_group]
  }
  beta_div_dist <- phyloseq::distance(physeq = physeq, method = method)
  beta_div_dist <- as(beta_div_dist, "matrix")
  dist_df <- data.frame()
  counter <- 1
  for (groups in names(group2samp)) {
    sub_dist <- beta_div_dist[group2samp[[groups]], group2samp[[groups]]]
    no_samp_col <- ncol(sub_dist)
    no_samp_row <- nrow(sub_dist)
    for (cols in seq(no_samp_col)) {
      if (cols > 1) {
        for (rows in seq((cols - 1))) {
          dist_df[counter, "sample_pair"] <- paste0(colnames(sub_dist)[cols], 
            "-", rownames(sub_dist)[rows])
          dist_df[counter, "group"] <- groups
          dist_df[counter, "beta_div_method"] <- method
          dist_df[counter, "beta_div_value"] <- sub_dist[rows, 
            cols]
          counter = counter + 1
        }
      }
    }
  }
  plot_boxplot <- ggplot(data = dist_df, aes(x = group, y = beta_div_value, 
    color = group)) + geom_boxplot(outlier.shape = NA) + 
    geom_jitter() + theme_bw() + xlab(group) + ylab(method) + 
    theme(axis.text.x = element_text(angle = 45, vjust = 1, 
      hjust = 1))
  list_Out <- list(data = dist_df, plot = plot_boxplot)
  return(list_Out)
}

Box and whisker plots of distance within group distances

#remotes::install_github("antonioggsousa/micrUBIfuns")
#library(micrUBIfuns)
T1_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Herbicide")
T1_beta_plot <- T1_beta$plot
T1_beta_plot <- T1_beta_plot + theme_classic()+ guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75)
T1_beta_plot

Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).

my_legend <- get_legend(T1_beta_plot)

Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).

as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_beta_legend.pdf")

Saving 7.29 x 4.51 in image

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_beta_legend.eps")

Saving 7.29 x 4.51 in image

T1_beta_plot<-T1_beta_plot+ theme(legend.position = "none") 
T1_beta_plot

Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T1_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image

Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T1_rare_withingroup_beta.eps")

Saving 7.29 x 4.51 in image

Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).

T1_beta_df<- T1_beta$data
T1_betamod<-aov(formula = beta_div_value ~ group ,data = T1_beta_df)
summary(T1_betamod)

             Df Sum Sq Mean Sq F value  Pr(>F)    
group         4 0.2967 0.07419    8.83 1.1e-06 ***
Residuals   250 2.1003 0.00840                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(x = T1_betamod, which = "group")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ group, data = T1_beta_df)

$group
                                        diff         lwr          upr     p adj
Dicamba-Atrazine-Mesotrione     -0.027993939 -0.07861966  0.022631781 0.5508308
Glyphosate-Atrazine-Mesotrione  -0.036412121 -0.08703784  0.014213599 0.2806138
Handweeded-Atrazine-Mesotrione  -0.058666667 -0.11176337 -0.005569963 0.0220819
Non-Treated-Atrazine-Mesotrione  0.040006061 -0.01061966  0.090631781 0.1940292
Glyphosate-Dicamba              -0.008418182 -0.05644596  0.039609594 0.9889807
Handweeded-Dicamba              -0.030672727 -0.08129845  0.019952993 0.4576855
Non-Treated-Dicamba              0.068000000  0.01997222  0.116027775 0.0012035
Handweeded-Glyphosate           -0.022254545 -0.07288027  0.028371175 0.7468226
Non-Treated-Glyphosate           0.076418182  0.02839041  0.124445957 0.0001747
Non-Treated-Handweeded           0.098672727  0.04804701  0.149298448 0.0000019

T2_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Herbicide")
T2_beta_plot <- T2_beta$plot
T2_beta_plot <- T2_beta_plot+ theme_classic() + theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75)
T2_beta_plot

Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 2 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T2_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image

Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 2 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T2_rare_withingroup_beta.eps")

Saving 7.29 x 4.51 in image

Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 2 rows containing missing values (`geom_point()`).

T2_beta_df<- T2_beta$data
T2_betamod<-aov(formula = beta_div_value ~ group ,data = T2_beta_df)
summary(T2_betamod)

             Df Sum Sq  Mean Sq F value  Pr(>F)   
group         4  0.080 0.019995   4.367 0.00201 **
Residuals   231  1.058 0.004578                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(x = T2_betamod, which = "group")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ group, data = T2_beta_df)

$group
                                        diff          lwr           upr     p adj
Dicamba-Atrazine-Mesotrione     -0.029188889 -0.070787897  0.0124101197 0.3048204
Glyphosate-Atrazine-Mesotrione  -0.020375758 -0.057770485  0.0170189703 0.5647456
Handweeded-Atrazine-Mesotrione   0.022206061 -0.015188667  0.0596007884 0.4780643
Non-Treated-Atrazine-Mesotrione -0.015555556 -0.054775477  0.0236643659 0.8113359
Glyphosate-Dicamba               0.008813131 -0.031069709  0.0486959716 0.9738007
Handweeded-Dicamba               0.051394949  0.011512109  0.0912777898 0.0043255
Non-Treated-Dicamba              0.013633333 -0.027965675  0.0552323419 0.8962631
Handweeded-Glyphosate            0.042581818  0.007106064  0.0780575719 0.0097763
Non-Treated-Glyphosate           0.004820202 -0.032574526  0.0422149299 0.9966025
Non-Treated-Handweeded          -0.037761616 -0.075156344 -0.0003668883 0.0465006

T3_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Herbicide") 
T3_beta$plot #+ scale_color_manual(values = c("#F8766D", "#A3A500",  "#00BF7D", "#00B0F6", "#E76BF3")) + 
T3_beta_plot <- T3_beta$plot
T3_beta_plot <- T3_beta_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("")
T3_beta_plot <-T3_beta_plot + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T3_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image

Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T3_rare_withingroup_beta.eps")

Saving 7.29 x 4.51 in image

Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).

T3_beta_df<- T3_beta$data
T3_betamod<-aov(formula = beta_div_value ~ group ,data = T3_beta_df)
summary(T3_betamod)

             Df Sum Sq  Mean Sq F value Pr(>F)
group         4 0.0672 0.016792   1.689  0.153
Residuals   252 2.5054 0.009942               

TukeyHSD(x = T3_betamod, which = "group")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ group, data = T3_beta_df)

$group
                                        diff         lwr        upr     p adj
Dicamba-Atrazine-Mesotrione      0.009452020 -0.04731207 0.06621612 0.9909320
Glyphosate-Atrazine-Mesotrione  -0.037518182 -0.08753732 0.01250096 0.2404830
Handweeded-Atrazine-Mesotrione  -0.014736364 -0.06475550 0.03528278 0.9275528
Non-Treated-Atrazine-Mesotrione -0.023409091 -0.07637301 0.02955483 0.7430333
Glyphosate-Dicamba              -0.046970202 -0.10570358 0.01176317 0.1840580
Handweeded-Dicamba              -0.024188384 -0.08292176 0.03454499 0.7896781
Non-Treated-Dicamba             -0.032861111 -0.09412180 0.02839957 0.5804264
Handweeded-Glyphosate            0.022781818 -0.02946147 0.07502511 0.7524585
Non-Treated-Glyphosate           0.014109091 -0.04096017 0.06917835 0.9554583
Non-Treated-Handweeded          -0.008672727 -0.06374199 0.04639653 0.9926689

library(ggpubr)
ggarrange(T1_beta_plot, T2_beta_plot, T3_beta_plot, ncol = 1)

Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_within_group_beta.pdf", width = 7, height = 10)

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_within_group_beta.eps", width = 7, height = 10)

Examination of dissimilarity across time points by treatment and then again by all chemical treatments combined.

T1_beta_df$Time<-"T1"
T2_beta_df$Time<-"T2"
T3_beta_df$Time<-"T3"


beta_div_T1_T2_T3 <- rbind(T1_beta_df, T2_beta_df, T3_beta_df)

beta_anova<-function(data, Herbicide = "Herbicide"){
  df_sub<- data %>% filter(group == Herbicide)
  mod<-aov(beta_div_value ~ Time, data = df_sub)
  print(summary(mod))
  print(TukeyHSD(mod, "Time"))
  boxplot(df_sub$beta_div_value ~ df_sub$Time)
}

#beta_anova(beta_div_T1_T2_T3, Herbicide = "Non-Treated")
#beta_anova(beta_div_T1_T2_T3, Herbicide = "Handweeded")
#beta_anova(beta_div_T1_T2_T3, Herbicide = "Dicamba")
#beta_anova(beta_div_T1_T2_T3, Herbicide = "Atrazine-Mesotrione")
#beta_anova(beta_div_T1_T2_T3, Herbicide = "Glyphosate")

#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide

index <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Handweeded", "Non-Treated")

sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])

#+ scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 


T1_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Mode")
T1_beta_chemical_combined_plot <- T1_beta_chemical_combined$plot 
T1_beta_chemical_combined_plot<- T1_beta_chemical_combined_plot + theme_classic() + guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 
T1_beta_chemical_combined_plot

Warning: Removed 9 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 9 rows containing missing values (`geom_point()`).

my_legend <- get_legend(T1_beta_chemical_combined_plot)

Warning: Removed 9 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 10 rows containing missing values (`geom_point()`).

as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_beta_combined_legend.pdf")

Saving 7.29 x 4.51 in image

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_beta_combined_legend.eps")

Saving 7.29 x 4.51 in image

T1_beta_chemical_combined_plot<-T1_beta_chemical_combined_plot+ theme(legend.position = "none")
T1_beta_chemical_combined_plot

Warning: Removed 9 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 10 rows containing missing values (`geom_point()`).

T2_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Mode")
T2_beta_chemical_combined_plot <- T2_beta_chemical_combined$plot 
T2_beta_chemical_combined_plot<- T2_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 
T2_beta_chemical_combined_plot

Warning: Removed 3 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 3 rows containing missing values (`geom_point()`).

T3_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Mode")
T3_beta_chemical_combined_plot <- T3_beta_chemical_combined$plot 
T3_beta_chemical_combined_plot<- T3_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 
T3_beta_chemical_combined_plot

Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).

ggarrange(T1_beta_chemical_combined_plot, T2_beta_chemical_combined_plot, T3_beta_chemical_combined_plot, ncol = 1)

Warning: Removed 9 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 9 rows containing missing values (`geom_point()`).
Warning: Removed 3 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 4 rows containing missing values (`geom_point()`).
Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_within_group_beta_chemical_combined.pdf", width = 5, height = 10)

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_within_group_beta_chemical_combined.eps", width = 5, height = 10)

T1_beta_df_chemical_combined <- T1_beta_chemical_combined$data
T2_beta_df_chemical_combined<- T2_beta_chemical_combined$data
T3_beta_df_chemical_combined<- T3_beta_chemical_combined$data

T1_beta_df_chemical_combined$Time<-"T1"
T2_beta_df_chemical_combined$Time<-"T2"
T3_beta_df_chemical_combined$Time<-"T3"

m1<-aov(beta_div_value ~ group, data = T1_beta_df_chemical_combined)
summary(m1)
TukeyHSD(m1, "group")
boxplot(beta_div_value ~ group, data = T1_beta_df_chemical_combined)


m2<-aov(beta_div_value ~ group, data = T2_beta_df_chemical_combined)
summary(m2)
TukeyHSD(m2, "group")
boxplot(beta_div_value ~ group, data = T2_beta_df_chemical_combined)

m3<-aov(beta_div_value ~ group, data = T3_beta_df_chemical_combined)
summary(m3)
TukeyHSD(m3, "group")
boxplot(beta_div_value ~ group, data = T3_beta_df_chemical_combined)


beta_div__chemical_combined_T1_T2_T3 <- rbind(T1_beta_df_chemical_combined, T2_beta_df_chemical_combined, T3_beta_df_chemical_combined)

beta_anova(beta_div__chemical_combined_T1_T2_T3, Herbicide = "Chemical")
beta_anova(beta_div__chemical_combined_T1_T2_T3, Herbicide = "Hand")
beta_anova(beta_div__chemical_combined_T1_T2_T3, Herbicide = "Non-Treated")

treatment to control

plotDistances = function(p, m, s, d) {

  # calc distances
  wu = phyloseq::distance(p, m)
  wu.m = melt(as.matrix(wu))
  
  # remove self-comparisons
  wu.m = wu.m %>%
    filter(as.character(Var1) != as.character(Var2)) %>%
    mutate_if(is.factor,as.character)
  
  # get sample data (S4 error OK and expected)
  sd = data.frame(sample_data(p)) %>%
    select(s, d) %>%
    mutate_if(is.factor,as.character)
  sd$Herbicide <- factor(sd$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
  
  # combined distances with sample data
  colnames(sd) = c("Var1", "Type1")
  wu.sd = left_join(wu.m, sd, by = "Var1")
  
  colnames(sd) = c("Var2", "Type2")
  wu.sd = left_join(wu.sd, sd, by = "Var2")
  
  #remove this line to plot all comparisons. 
  wu.sd = wu.sd %>% filter(Type1 == "Hand" | Type1 == "Non-Treated")
  
  # plot
  ggplot(wu.sd, aes(x = Type2, y = value)) +
    theme_bw() +
    geom_point() +
    geom_boxplot(aes(color = ifelse(Type1 == Type2, "red", "black"))) +
    scale_color_identity() +
    facet_wrap(~ Type1, scales = "free_x") +
    theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
    ggtitle(paste0("Distance Metric = ", m))
  
}

a<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T1"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
a <- a + ggtitle("Time 1 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T1_rare_allgroup_beta.pdf")
b<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T2"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
b <-b + ggtitle("Time 2 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T2_rare_allgroup_beta.pdf")
c<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T3"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
c<- c + ggtitle("Time 3 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T3_rare_allgroup_beta.pdf")

library(ggpubr)
ggarrange(a, b, c, ncol = 1)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_allgroup_beta.pdf", width = 7, height = 12)

Taxon abundance bar plot

#create super long color vector
col_vector <- c("#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941", "#006FA6", "#A30059",
        "#FFDBE5", "#7A4900", "#0000A6", "#63FFAC", "#B79762", "#004D43", "#8FB0FF", "#997D87",
        "#5A0007", "#809693", "#FEFFE6", "#1B4400", "#4FC601", "#3B5DFF", "#4A3B53", "#FF2F80",
        "#61615A", "#BA0900", "#6B7900", "#00C2A0", "#FFAA92", "#FF90C9", "#B903AA", "#D16100",
        "#DDEFFF", "#000035", "#7B4F4B", "#A1C299", "#300018", "#0AA6D8", "#013349", "#00846F",
        "#372101", "#FFB500", "#C2FFED", "#A079BF", "#CC0744", "#C0B9B2", "#C2FF99", "#001E09",
        "#00489C", "#6F0062", "#0CBD66", "#EEC3FF", "#456D75", "#B77B68", "#7A87A1", "#788D66",
        "#885578", "#FAD09F", "#FF8A9A", "#D157A0", "#BEC459", "#456648", "#0086ED", "#886F4C",
        
        "#34362D", "#B4A8BD", "#00A6AA", "#452C2C", "#636375", "#A3C8C9", "#FF913F", "#938A81",
        "#575329", "#00FECF", "#B05B6F", "#8CD0FF", "#3B9700", "#04F757", "#C8A1A1", "#1E6E00",
        "#7900D7", "#A77500", "#6367A9", "#A05837", "#6B002C", "#772600", "#D790FF", "#9B9700",
        "#549E79", "#FFF69F", "#201625", "#72418F", "#BC23FF", "#99ADC0", "#3A2465", "#922329",
        "#5B4534", "#FDE8DC", "#404E55", "#0089A3", "#CB7E98", "#A4E804", "#324E72", "#6A3A4C",
        "#83AB58", "#001C1E", "#D1F7CE", "#004B28", "#C8D0F6", "#A3A489", "#806C66", "#222800",
        "#BF5650", "#E83000", "#66796D", "#DA007C", "#FF1A59", "#8ADBB4", "#1E0200", "#5B4E51",
        "#C895C5", "#320033", "#FF6832", "#66E1D3", "#CFCDAC", "#D0AC94", "#7ED379", "#012C58",
        
        "#7A7BFF", "#D68E01", "#353339", "#78AFA1", "#FEB2C6", "#75797C", "#837393", "#943A4D",
        "#B5F4FF", "#D2DCD5", "#9556BD", "#6A714A", "#001325", "#02525F", "#0AA3F7", "#E98176",
        "#DBD5DD", "#5EBCD1", "#3D4F44", "#7E6405", "#02684E", "#962B75", "#8D8546", "#9695C5",
        "#E773CE", "#D86A78", "#3E89BE", "#CA834E", "#518A87", "#5B113C", "#55813B", "#E704C4",
        "#00005F", "#A97399", "#4B8160", "#59738A", "#FF5DA7", "#F7C9BF", "#643127", "#513A01",
        "#6B94AA", "#51A058", "#A45B02", "#1D1702", "#E20027", "#E7AB63", "#4C6001", "#9C6966",
        "#64547B", "#97979E", "#006A66", "#391406", "#F4D749", "#0045D2", "#006C31", "#DDB6D0",
        "#7C6571", "#9FB2A4", "#00D891", "#15A08A", "#BC65E9", "#FFFFFE", "#C6DC99", "#203B3C",

        "#671190", "#6B3A64", "#F5E1FF", "#FFA0F2", "#CCAA35", "#374527", "#8BB400", "#797868",
        "#C6005A", "#3B000A", "#C86240", "#29607C", "#402334", "#7D5A44", "#CCB87C", "#B88183",
        "#AA5199", "#B5D6C3", "#A38469", "#9F94F0", "#A74571", "#B894A6", "#71BB8C", "#00B433",
        "#789EC9", "#6D80BA", "#953F00", "#5EFF03", "#E4FFFC", "#1BE177", "#BCB1E5", "#76912F",
        "#003109", "#0060CD", "#D20096", "#895563", "#29201D", "#5B3213", "#A76F42", "#89412E",
        "#1A3A2A", "#494B5A", "#A88C85", "#F4ABAA", "#A3F3AB", "#00C6C8", "#EA8B66", "#958A9F",
        "#BDC9D2", "#9FA064", "#BE4700", "#658188", "#83A485", "#453C23", "#47675D", "#3A3F00",
        "#061203", "#DFFB71", "#868E7E", "#98D058", "#6C8F7D", "#D7BFC2", "#3C3E6E", "#D83D66",

        "#2F5D9B", "#6C5E46", "#D25B88", "#5B656C", "#00B57F", "#545C46", "#866097", "#365D25",
        "#252F99", "#00CCFF", "#674E60", "#FC009C", "#92896B")

phylumGlommed <- tax_glom(ps_rare, "Phylum")

#t1
phylumGlommed_herb_t1 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T1"), group = "Herbicide")
phylumGlommed_herb_t1 <- transform_sample_counts(phylumGlommed_herb_t1, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t1)$Herbicide <- factor(sample_data(phylumGlommed_herb_t1)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t1, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_Taxon_barplot_t1.pdf")

#t2
phylumGlommed_herb_t2 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T2"), group = "Herbicide")
phylumGlommed_herb_t2 <- transform_sample_counts(phylumGlommed_herb_t2, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t2)$Herbicide <- factor(sample_data(phylumGlommed_herb_t2)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t2, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/ITS_Taxon_barplot_t2.pdf")

#t3
phylumGlommed_herb_t3 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T3"), group = "Herbicide")
phylumGlommed_herb_t3 <- transform_sample_counts(phylumGlommed_herb_t3, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t3)$Herbicide <- factor(sample_data(phylumGlommed_herb_t3)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t3, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/ITS_Taxon_barplot_t3.pdf")

Combined herbicide and time bar plot

sample_data(ps_rare_sub)$herb_time<-paste(sample_data(ps_rare_sub)$Herbicide, sample_data(ps_rare_sub)$Time, sep = "_")
ps_rare_for_barplot <- prune_taxa(taxa_sums(ps_rare_sub) > 50, ps_rare_sub)
plot_bar(ps_rare_for_barplot, x= "herb_time", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_BarPlot_Herbicide_Time.pdf", width = 20, height = 11)

Linear modeling of abundant taxa

mod_abund(test, IV = "Time")

                   Df Sum Sq Mean Sq F value Pr(>F)  
matrix(data[, IV])  2  30.58  15.292   3.725 0.0378 *
Residuals          26 106.73   4.105                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
[1] "Dicamba"           "f.Coniochaetaceae"
                   Df Sum Sq Mean Sq F value  Pr(>F)   
matrix(data[, IV])  2  864.4   432.2   7.708 0.00207 **
Residuals          29 1626.0    56.1                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
[1] "Handweeded"           "f.Spizellomycetaceae"

                   Df Sum Sq Mean Sq F value Pr(>F)  
matrix(data[, IV])  2    805   402.7   3.351 0.0496 *
Residuals          28   3366   120.2                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
[1] "Non-Treated"          "f.Spizellomycetaceae"

                   Df Sum Sq Mean Sq F value Pr(>F)  
matrix(data[, IV])  2  397.2  198.58   4.159 0.0255 *
Residuals          30 1432.4   47.75                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
[1] "Glyphosate"          "f.Ceratobasidiaceae"

                   Df Sum Sq Mean Sq F value   Pr(>F)    
matrix(data[, IV])  2  204.1  102.03   9.469 0.000649 ***
Residuals          30  323.3   10.78                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
[1] "Glyphosate"           "f.Spizellomycetaceae"